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Abstract — The  field  of  multi-target  tracking  faces  two  primary 
challenges:  (i)  data  association  and  (ii)  trajectory  estimation. 
MTT  problems  are  well  researched  with  many  algorithms  solving 
these  two  problems  separately,  however  few  algorithms  attempt  to 
solve  these  simultaneously  and  even  fewer  utilize  optimization.  In 
this  paper  we  introduce  a  new  mixed  integer  optimization  (MIO) 
model  which  solves  the  data  association  and  trajectory  estimation 
problems  simultaneously  by  minimizing  an  easily  interpretable 
global  objective  function.  Furthermore,  we  propose  a  greedy 
heuristic  which  quickly  finds  good  solutions.  We  extend  both  the 
heuristic  and  the  MIO  model  to  scenarios  with  missed  detections 
and  false  alarms. 

Index  Terms — optimization;  multi-target  tracking;  data  asso¬ 
ciation;  trajectory  estimation;  mixed  integer  optimization 

I.  Introduction 

ULTI-target  tracking  is  the  problem  of  estimation  the 
state  of  multiple  dynamic  objects,  referred  to  as  targets 
over  a  fixed  window  of  time.  At  various  points  of  time  within 
the  window,  the  targets  are  observed  in  a  scan,  resulting  in  set 
of  detections.  From  these  detections,  the  multi-target  tracking 
problem  aims  to  extract  information  about  target  dynamics. 

Solutions  to  this  problem  are  sought  across  many  civilian 
and  military  applications  including  but  not  limited  to  ballistic 
missile  and  aircraft  defense,  space  applications,  the  movement 
of  ships  and  ground  troops,  autonomous  vehicles  and  robotics, 
and  air  traffic  control.  Each  application  has  unique  attributes 
and  assumptions,  and  various  algorithms  have  been  developed 
for  each.  As  a  result,  the  field  of  multi-target  tracking  has 
expanded  to  numerous  research  venues,  and  there  is  a  wide 
range  of  literature  on  the  topic.  A  more  complete  overview 
of  all  MTT  methods,  including  the  classes  of  algorithms  and 
their  variants  as  well  as  additional  methods  not  discussed  in 
this  paper,  can  be  found  in  [1],  For  a  more  exhaustive  overview 
of  estimation  techniques,  filtering,  gating,  and  more  please  see 
[2]  or  [3], 

The  field  of  multi-target  tracking  faces  two  primary  chal¬ 
lenges:  (i)  data  association  and  (ii)  trajectory  estimation.  Given 
a  set  of  sensor  detections  the  data  association  problem  consists 
of  assigning  the  detections  to  a  set  of  targets.  Alternatively, 
this  can  be  viewed  as  a  labeling  problem  in  which  each 
detection  needs  to  be  labeled  with  a  target  identifier.  The 
association  problem  is  further  complicated  when  sensors  fail 
to  report  detections  (missed  detection)  or  incorrectly  report 
detections  (false  alarm),  resulting  in  ambiguity  in  the  number 
of  existing  targets.  The  trajectory  estimation  problem  consists 
of  estimating  the  state  space  of  a  target  (i.e.,  position,  velocity, 
acceleration,  size,  etc.)  from  the  associated  detections  of  the 
aforementioned  assignment  problem.  Even  when  all  of  the 


associations  are  known,  the  estimation  problem  is  challenging 
due  to  the  presence  of  measurement  noise.  As  can  be  seen, 
the  two  problems  of  data  association  and  trajectory  estimation 
are  closely  related  and  dependent  on  one  another. 

Some  classical  algorithms  treat  the  data  association  and 
trajectory  estimation  problems  separately  using  a  combination 
of  probabilistic  approaches  to  determine  data  associations  and 
filters  to  estimate  trajectories.  One  such  algorithm  is  the  global 
nearest  neighbor  (GNN).  The  GNN  algorithm  is  a  naive  2-D 
assignment  algorithm,  which  evaluates  one  scan  of  detections 
at  a  time,  globally  assigning  the  nearest  detection  at  each 
scan  [4],  Once  the  data  association  has  been  determined, 
the  detections  are  often  passed  through  one  of  numerous 
filters,  most  commonly  a  Kalman  filter  [5],  which  updates  the 
trajectory  estimates  before  the  algorithm  progresses  forward  to 
the  next  scan.  This  process  repeats  sequentially  through  each 
scan  of  data. 

Modern  algorithms  in  the  field  of  multi-target  tracking 
are  most  commonly  statistical  based,  often  relying  on  heavy 
probabilistic  assumptions  about  the  underlying  target  dynam¬ 
ics  or  detection  process.  The  two  most  prevalent  statistical 
algorithms  in  the  field  of  multi-target  tracking  are  the  the 
Multiple  Hypothesis  Tracker  (MHT)  and  the  Joint  Probability 
Data  Association  Filter  (JPDAF)  and  their  numerous  variants 
and  extensions.  Both  classes  of  algorithms  attempt  to  solve 
the  data  association  problem  by  generating  a  set  of  potential 
hypotheses,  or  possible  detection-to-track  assignments.  Here  a 
track  is  a  set  of  labelled  detections  belonging  to  the  same 
target.  Probabilities  are  assigned  to  each  hypothesis  based 
on  the  likelihood  of  the  trajectory’s  existence,  and  numerous 
approaches  for  accomplishing  this  task  have  been  proposed. 

The  MHT,  first  proposed  by  Reid  in  [6],  assigns  likelihood 
values  to  hypotheses  using  a  Bayesian  MAP  estimator,  which 
requires  assumptions  on  object  dynamics.  This  algorithm  is 
generally  considered  to  be  the  modern  standard  for  solving 
the  data  association  problem.  Many  variants  have  been  pro¬ 
posed  for  implementation  which  leverage  techniques  such  as 
clustering,  gating,  hypothesis  selection,  hypothesis  pruning, 
and  merging  of  state  estimates.  Many  of  these  methods  are 
summarized  by  Blackman  in  [7], 

While  the  MHT  has  seen  various  forms  of  success,  it  faces 
several  key  challenges.  Namely,  the  curse  of  dimensionality 
and  complexity.  The  number  of  possible  hypotheses  grows 
exponentially  with  the  number  of  potential  tracks  and  the 
number  of  scans.  Consequently,  it  is  considered  intractable  for 
large  scenarios.  Moreover,  the  MHT  might  require  extensive 
tuning  and  thus  may  be  difficult  to  implement  in  practice,  in 
addition  to  being  computationally  expensive.  For  these  reasons 


2 


it  is  generally  considered  to  be  one  of  the  most  complex  MTT 
algorithms. 

A  Probability  Data  Association  (PDA)  takes  a  Bayesian 
approach  to  solving  the  data  association  problem  by  finding 
detection-to-target  assignment  probabilities  via  a  posterior 
PDF,  which  again  requires  heavy  assumptions  on  object  dy¬ 
namics  and  the  detection  process.  In  similar  fashion,  a  Joint 
PDA  (JPDA)  assigns  probabilities  that  are  computed  jointly 
across  all  targets.  The  JPDAF  is  an  algorithm  which  imple¬ 
ments  the  JPDA  along  with  filters  and  estimation  methods  as 
discussed  previously.  [2] 

A  limited  number  of  optimization  based  algorithms  have 
been  applied  to  solve  the  MTT  problem,  most  of  which  attempt 
to  solve  by  mapping  the  measurement  set  onto  a  trellis  and 
seek  the  optimal  measurement  association  sequence.  Some 
examples  include  the  Multi-Target  Viterbi[8]  and  an  extension 
in  [9]  which  formulates  [8]  as  a  network  flow,  reducing 
the  solve  time  from  exponential  to  polynomial.  Others  have 
suggested  adaptations  which  allow  this  approach  to  be  used 
similar  to  MHT  methods  by  outputting  a  single  best  set  of  K 
tracks,  or  a  list  of  L  best  sets  of  k  tracks  [10]. 

Compared  to  the  number  of  statistical  based  algorithms 
in  the  MTT  literature,  optimization  based  algorithms  are 
relatively  lacking.  In  fact,  most  occurrences  of  optimization 
in  the  MTT  literature  propose  the  use  of  optimization  to 
leverage  statistical  algorithms,  in  particular  the  MHT.  For 
example.  Integer  optimization  has  been  used  to  improve  MHT 
hypothesis  selection  by  solving  an  assignment  problem  which 
chooses  the  best  hypothesis,  but  only  after  costs  have  been 
assigned  (statistically  based)  and  hypotheses  have  been  pruned 

[11] ,  Somewhat  similarly,  linear  optimization  has  also  been 
used  to  assist  in  the  hypothesis  selection  process  for  the  MHT 

[12] ,  Still,  other  attempts  aim  to  improve  the  MHT  hypothesis 
selection  process  via  Lagrangian  relaxation  [13]. 

More  recently,  Andriyenko  and  Schindler  have  proposed 
formulating  the  MTT  problem  as  a  minimization  of  a  con¬ 
tinuous  energy  in  [14]  and  then  again  as  a  minimization 
of  discrete-continuous  energy  in  [15],  These  algorithms  aim 
to  more  accurately  represent  the  nature  of  the  problem,  but 
sacrifice  interpretability  for  complexity  in  the  process.  Rather 
than  formulating  the  problem  to  lend  it  easily  to  traditional 
global  optimization  methods,  the  authors  intend  to  leverage  the 
use  of  optimization  techniques  to  find  strong  local  minima  of 
their  proposed  energy  objective,  and  they  achieve  strong  results 
in  doing  so.  However,  this  approach  calls  for  the  use  of  several 
parameters  that  must  be  tuned  and  few  recommendations 
are  provided  for  how  to  go  about  such  a  tuning  process. 
Additionally,  these  methods  require  initialization  heuristics  to 
begin  the  solving  process,  which  is  in  itself  complicated  to 
implement  and  is  not  directly  connected  to  the  optimization 
problem  solved. 

In  this  paper  we  propose  the  use  of  mixed  integer  op¬ 
timization  (MIO)  to  formulate  and  solve  the  multi-target 
tracking  problem.  Although  MIOs  are  generally  thought  to  be 
intractable  (NP-Hard),  in  many  practical  cases  near  optimal 
solutions  and  even  optimal  solutions  to  these  problems  can 
be  obtained  in  reasonable  time  [16].  This  can  be  attributed  to 
the  fact  that  MIO  solvers  have  seen  significant  performance 


improvements  in  recent  years  due  to  advancements  in  both 
methodology  and  hardware.  The  development  of  new  heuristic 
methods,  discoveries  in  cutting  plane  theory,  and  improved 
linear  optimization  methods  have  all  contributed  to  improve¬ 
ments  in  performance  [17].  Modern  solvers  such  as  Gurobi 
and  CPLEX  have  been  shown  to  perform  extremely  well  on 
benchmark  tests.  In  the  past  six  years  alone,  Gurobi  has  seen 
performance  improvements  by  a  factor  of  48.7  [18].  CPLEX 
saw  improvements  by  a  factor  of  29,000  from  1991  to  2007 
[19],  From  1994  to  2014,  the  growth  of  supercomputing  power 
as  recorded  by  the  TOP500  list  has  improved  by  a  factor  of 
567839  [20].  Thus,  the  total  combined  effective  improvement 
of  software  and  hardware  advancements  is  on  the  scale  of  800 
billion  times  in  the  past  25  years. 

The  literature  is  also  lacking  in  performance  metrics  for  the 
evaluation  of  MTT  algorithms.  There  is  no  standard  method  of 
measuring  scenario  complexity  or  algorithm  performance  as  a 
function  of  this  complexity.  In  many  cases  only  the  sensor’s 
detection  noise  is  taken  into  account  and  other  factors  such 
as  target  density  is  negated.  Recent  work  [21]  proposes  a 
mathematically  rigorous  performance  metric  for  measuring  the 
distance  between  ground  truth  and  estimated  track,  but  there  is 
not  much  attention  given  to  the  complexity  of  generated  sce¬ 
narios.  In  this  paper  we  also  suggest  measures  of  complexity 
and  performance  which  are  related  to  the  ones  suggested  in 
[2 1  ]  but  we  show  the  value  in  relating  a  complexity  measure  to 
performance  measures,  namely  that  it  allows  you  to  evaluate 
the  data  association  and  trajectory  estimation  problems  sepa¬ 
rately.  We  evaluate  the  methods  suggested  in  this  paper  using 
these  complexity  and  performance  measures  on  two  simulated 
experiments. 

The  main  contributions  of  this  paper  are  as  follows: 

(i)  We  introduce  a  simple  interpretable  MIO  model  which 
solves  the  data  association  and  trajectory  estimation 
problems  simultaneously  for  a  sensor  with  no  detection 
ambiguity.  The  model  does  not  require  to  tuning  of 
parameters.  This  MIO  is  tractable,  in  the  sense  that  it  can 
be  solved  to  optimality  or  near  optimality  in  a  reasonable 
amount  of  time,  for  the  considered  applications. 

(ii)  We  propose  a  heuristic,  motivated  by  the  optimization 
problem,  which  gives  us  feasible  solutions  to  this  prob¬ 
lem  and  show  how  it  can  be  used  as  warm  start  to  the 
MIO  in  order  to  improve  the  quality  of  the  solutions 
obtained  as  well  as  the  running  time. 

(iii)  We  extend  this  basic  MIO  model  and  corresponding 
heuristic  initialization  algorithm  for  the  case  of  detection 
ambiguity,  i.e.,  the  case  where  there  are  both  missed  de¬ 
tections  and  false  alarms,  keeping  interpretability  while 
only  adding  two  tunable  parameters,  as  well  as  provide 
general  guidelines  as  to  how  tune  this  parameters. 

(iv)  We  present  a  new  measure  of  complexity  for  the  data 
association  problem,  and  show  how  it  aids  in  scenario 
generation.  We  also  discuss  a  simplified  measure  of 
performance  for  the  trajectory  estimation  problem. 

The  paper  structure  is  as  follows.  We  begin  with  a  de¬ 
scription  of  the  MTT  problem  as  we  wish  to  model  it  in 
Section  II.  In  Section  III  we  develop  a  simple  MIO  formulation 
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for  a  sensor  with  no  detection  ambiguity  and  extend  it  to 
a  generalized  formulation.  Following  is  a  discussion  on  a 
proposed  heuristic  in  Section  IV.  Then  we  propose  extensions 
for  both  the  MIO  formulation  and  the  heuristic  to  a  sensor 
with  detection  ambiguity  in  Section  V.  Metrics  for  measuring 
scenario  complexity  and  algorithm  performance  are  proposed 
in  Section  VI.  Experimental  simulations  are  outlined  in  Sec¬ 
tion  VII.  A  summary  of  significant  computational  results  are 
discussed  in  Section  VIII.  Finally,  conclusions  and  future  work 
are  discussed  in  Section  IX. 

General  Notations:  Unless  specified  otherwise,  ||  •  |  is  used 
to  indicate  the  LI  norm,  and  |  •  |  refers  to  element  wise  absolute 
value. 


Notation:  We  observe  P  targets  over  a  fixed  time  window 
over  which  T  scans  are  collected.  Scans  occur  at  a  fixed  rate, 
usually  of  about  1Hz,  such  that  the  set  of  scans  is  denoted  by 
{ti,  t-2,  The  ith  detection  of  the  tth  scan  is  indicated  by 

Xu,  such  that  a  scan  of  data  at  time  t  is  the  unordered  set  of 
detections  Xt  =  {in,  X2,t,  •••>  xp, t}-  The  data  for  the  problem 
is  the  ordered  set  of  scans  X  =  {X±,  X2, Xp}.  The  state 
space  of  target  trajectories  is  paramatarized  by  a  true  initial 
position  a‘rue  and  a  true  constant  velocity  /3jr“e.  Therefore, 
the  true  position  Xjt  of  trajectory  j  at  scan  t  is  given  by: 

xjt  =  afue  +  pfuet  (1) 


II.  Problem  Description 

In  this  paper,  we  restrict  our  exploration  of  the  MTT 
problem  to  the  automatic  tracking  of  multiple,  independent 
point  targets  using  a  single  sensor.  A  target  is  the  object  of 
interest.  A  point  target’s  only  identifiable  attributes  are  its  state 
space,  which  we  restrict  to  position  and  velocity.  The  state 
space  fully  defines  the  field  of  trajectories,  or  paths  along 
which  targets  travel.  A  detection  is  collected  from  each  target 
at  sequential  scans.  Detections  are  subject  to  noise.  We  treat 
two  general  scenarios:  with  and  without  detection  ambiguity. 

When  there  is  no  detection  ambiguity,  the  sensor  produces 
exactly  one  detection  for  each  target  at  each  time,  and  there 
is  no  other  source  of  detections.  Therefore,  the  number  of 
detections  at  each  point  in  time  is  the  same  as  the  number 
of  targets,  and  the  data  association  problem  at  each  point  in 
time  is  equivalent  to  a  simple  assignment  problem.  Our  basic 
optimization  model,  presented  in  section  [add  reference]  will 
this  with  this  problem. 

Detection  ambiguity  refers  to  the  more  complex  case  where 
the  sensor  generates  both  false  alarms  and  missed  detections. 
A  False  Alarm  occurs  when  a  detection  is  collected  when  no 
target  exists.  This  could  be  the  result  of  measurement  error 
or  difficulties  in  signal  processing.  A  Missed  Detection  occurs 
when  a  data  point  is  not  collected  at  a  given  time  when  a 
target  actually  exists.  Therefore,  the  number  of  detections  at 
each  point  in  time  may  be  either  higher  or  lower  than  the 
actual  number  of  targets,  and  each  detection  can  be  classified 
in  either  of  these  categories  in  addition  to  assigning  targets 
to  trajectories  as  before.  In  section  [add  reference]  we  will 
extend  the  formulation  of  basic  model  to  a  robust  formulation 
dealing  with  this  ambiguity,  and  refer  to  it  as  the  robust  MIO 
model. 

Throughout  the  paper  we  make  the  following  assumptions: 

Assumption  1.  (i)  All  targets  have  constant  velocity,  i.e.. 

Targets  do  not  maneuver  and  no  outside  forces  act  on 
them. 

(ii)  Each  target’s  dynamics  are  independent  of  one  another. 
(Hi)  The  number  of  targets  remains  constant  throughout  the 
window  of  observation,  i.e.,  there  is  no  birth/death  of 
targets. 

(iv)  Each  target  produces  at  most  one  detection  per  scan. 

(v)  The  detection  errors  are  independent  of  one  another. 


III.  Basic  MIO  Model 

In  this  section,  we  deal  with  the  case  of  no  detection 
ambiguity.  Therefore,  we  add  the  following,  more  restrictive 
assumptions,  to  those  presented  in  Assumption  1 

Assumption  2.  (i)  The  sensor  generates  exactly  one  detec¬ 

tion  for  each  target  at  each  time  (no  missed  detections), 
(ii)  The  sensor  does  not  generate  any  additional  detections 
(no  false  alarms). 

We  begin  constructing  our  MIO  model  by  defining  decision 
variables  that  represent  the  desired  detection  to  target  asso¬ 
ciations  and  target  estimated  trajectories.  Next,  using  these 
decision  variables,  we  develop  an  objective  function  which 
mathematically  quantifies  the  value  of  the  model  decisions,  in 
this  case  as  a  measure  of  distance  of  the  estimated  trajectories 
from  the  associated  detections.  Finally,  we  restrict  these  vari¬ 
ables  using  constraints  that  force  the  model  to  find  solutions 
that  are  feasible  for  the  MTT  problem.  The  model  is  developed 
step  by  step  in  the  coming  sections  before  the  full  model  is 
presented. 


A.  Decision  Variables 

The  data  association  and  trajectory  estimation  problems 
require  unique  decision  variables.  Because  these  two  problems 
lie  in  different  domains,  the  variables  we  use  to  represent 
these  decisions  also  differ.  First,  we  introduce  continuous 
decision  variables  ay  £  Rn  and  f3j  £  R"  to  represent  the 
estimated  initial  position  and  velocity  of  each  trajectory  j.  In 
our  interpretation  of  the  MTT  problem  we  allow  the  trajectory 
parameters  to  lie  anywhere  in  the  real-continuous  domain.  For 
the  data  estimation  problem,  we  wish  to  assign  detections 
to  trajectories,  a  naturally  discrete  problem.  Therefore,  we 
introduce  binary  decision  variables  yitj  to  indicate  whether 
detection  x,f  is  assigned  to  trajectory  j  or  not: 


if  detection  Xu  is  assigned  to  trajectory  j 
otherwise. 


’  (2) 


B.  Objective  Function 

Next,  we  would  like  to  develop  a  function  which  accurately 
scores  the  quality  of  a  feasible  solution.  An  ideal  objective 
function  would  jointly  provide  a  single  quantifiable  measure 
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of  goodness  for  both  the  data  association  and  trajectory 
estimation  problems.  Therefore  we  want  the  objective  to  take 
into  account  the  assignments  of  detections  in  addition  to  the 
estimated  trajectory  determined  by  those  assignments. 

In  terms  of  our  established  decision  variables,  the  estimated 
position  of  the  linear  trajectory  j  at  time  t  is  given  by: 

Xjt  =  dj  +  /3jt  (3) 


C.  Constraints 

For  each  scan,  each  detection  Xu  must  be  assigned  to 
exactly  one  target  j: 

p 

E  Uitj  =  1  Vi,  t  (10) 

j=i 

Similarly,  for  each  scan,  each  target  must  be  assigned 
exactly  one  detection: 


For  the  trajectory  estimation  problem  we  wish  to  minimize 
the  distance  between  x,;4  and  Xjt.  In  other  words,  for  detections 
Xu  assigned  to  trajectory  j,  we  wish  to  minimize  \\xu  —  Xjt\\ 
for  some  norm.  The  two  natural  norms  to  consider  here  are  the 
LI  and  L2  norms.  The  LI  norm  has  the  advantage  that  it  can  be 
reformulated  using  linear  optimization  (through  the  addition  of 
continuous  variables  and  constraints),  and  it  is  well  known  to 
be  more  robust  to  outliers.  Additionally,  existing  algorithms  for 
MIO  are  more  well  developed  for  linear  rather  than  quadratic 
optimization.  However,  the  L2  norm  square  form  (RSS)  has 
the  advantage  that  it  can  be  quickly  computed  using  a  matrix 
formulation,  making  it  more  predisposed  to  a  heuristic.  This 
concept  will  be  discussed  further  in  section  IV. 

Substituting  (3)  for  Xjt  we  arrive  at  our  objective  function: 

T 

minimize:  E  E  II Xit  -  dj  -  (3jt\\  (4) 

where  A  is  the  set  of  pairs  that  indicate  the  assignment  of 
detection  i  to  trajectory  j. 

For  the  data  association  problem,  we  wish  to  only  penalize 
the  objective  function  when  detection  xit  has  been  assigned 
to  trajectory  j,  which  occurs  when  ynj  =  1.  An  easy  method 
to  enforce  this  using  our  established  variables  would  be  to 
construct  an  interaction  term  like  in  (5)  below. 

p  T 

minimize:  Y^  Y^  \yujXit  -  dj  -  /3jt\\  (5) 

We  now  show  that  this  objective  can  be  converted  to  a  linear 
program  in  the  case  of  the  LI  norm  by  introducing  continuous 
variables  9jt  and  the  following  additional  constraints. 


p 

'S^Vitj  =  1  V?,f  (11) 

i= 1 


D.  Simple  Formulation 

Combining  all  of  these  elements  together,  we  arrive  at  the 
following  MIO  model: 


minimize: 

Ojt 


subject  to: 


p  T 

E  E 

3= 1  t=l 
P 

'y '  yaj  —  i  f 
j= i 

p 

y  ' yaj =  i  yjit 

i—1 

yaj  Xu  dj  (3jt  G  0 jt 

{VitjXit  dj  ly  t )  L  9 jt  VL  j1" ,  t 

yitj&{  0,1}  Vi,t,3 

dj  G  M"  Vj,  pj  G  Mn  Vj,  Zjt  G  R" 


Vj,f 


E.  Generalized  Formulation 

The  previous  formulation,  although  simple  and  easily  inter¬ 
pretable,  has  the  disadvantages  of  being  (i)  dense  and  (//)  ill 
suited  for  extension  to  detection  ambiguity.  Therefore,  next, 
we  present  a  generalized  formulation  which  can  be  linearized 
through  the  introduction  of  additional  binary  decision  vari¬ 
ables.  Alternatively,  we  can  create  a  new  variable  Zjt  which 
takes  on  the  value  Xu  when  ijijt  =  1  and  some  arbitrary 
number  when  ynj  =  0.  Using  this  method  we  must  adjust 
the  objective  function  below. 


VitjXit  dj  (3jt  G  9 jt  Vi,  j,  t  (6) 

( VitjXit  dj  l3j  t )  G  Ojt  V/" ,  j .  /  (7) 

The  resulting  objective  function  for  the  case  of  the  LI  norm 
would  then  be: 


p  T 

linimize:  E  E  On 

3t  j= i  t=i 


minimize 


p  T 

■  EE 

i= 1 i=l 


9jt  ||  2 


P  T 


minimize:  2^  E  “  ai  ~  PM  (12) 

j- 1 1= i 

This  objective  can  then  be  linearized  by  again  introducing 
Ojt  and  similar  constraints  as  follows. 


(8) 


where  e  is  the  vector  of  ones. 

For  the  case  of  the  L2  norm,  the  objective  function  would 
be: 


p  T 

minimize: 

3=1  i= 1 


e 


jt 


Zj t  dj  3 j i  G  o, t  Vi,j,t 

dj  3j  t )  3  Ojt  Vi,j,  t 


(13) 

(14) 

(15) 


(9) 


Furthermore,  we  must  ensure  that  the  decision  variable  Zjt 
will  only  take  on  the  value  of  xif  in  the  objective  function  if 
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Xu  is  assigned  to  target  j  (yuj  =  1).  We  enforce  this  effect 
using  the  following  constraint: 

HR ( 1  Vitj)  —  | Zjt  |  Vt,f,  j  (16) 

(17) 

where  Mt  =  max|a;jt|  for  each  scan.  Furthermore,  we  can 

j 

write  this  equivalently  as  a  linear  optimization  problem  by 
using  the  following  set  of  two  linear  constraints: 

'EitVitj  T  (  1  Vitj )  ^  Zjt  V? ■  t ■  j  (18) 

% itUitj  HRf 1  Vitj')  V  7y /,  Vi •  t • j  (19) 

Combining  all  of  these  elements  together,  we  arrive  at  the 
following  generalized  MIO  model: 

P  T 

minimize:  E  E 

j=i  t=i 
P 

subject  to:  E  Vitj  1  Vi ,  t 

3= 1 
P 

E  yuj  =  i  vj,  t 

i= 1 

itUitj  '  A7t(  1  Vitj')  V  Vi ,  t  ■  J 

^itUitj  V/f  ( 1  Vitj)  —  Zjt  Vi  ,  I,  j 

Zjt  rty  {3 j  1.  G  0:l f  Vi,  j,  t 

(Vp  uy  V/ 1  )  ^  0 jt  Vi,  j,  t 

Vitj  &  {0,1}  Vi,t,j 

otj  G  Rn  Vj,  ^  G  Rn  Vj,  G  Rn  Vj,  f 

IV.  Heuristic 

Next,  we  present  a  detailed  description  of  a  heuristic  which 
finds  good  feasible  solutions.  These  solutions  can  be  used  as  a 
warm  start  to  the  MIO,  providing  a  performance  boost  to  the 
MIO.  The  heuristic  leverages  the  power  of  randomization  local 
search  methods  to  find  locally  optimal  solutions.  Although 
the  heuristic  takes  a  local  search  approach,  we  hypothesize 
that  it  will  discover  near  optimal  solutions,  provide  that  it  is 
initialized  with  enough  random  starting  points. 

An  important  distinction  to  discuss  here  is  the  difference 
in  objective  functions.  As  discussed  before  the  two  natural 
choices  are  the  LI  and  L2  norms.  For  the  heuristic,  we  desire 
an  objective  which  can  be  calculated  efficiently.  Therefore, 
in  this  case  the  L2  norm  square  (RSS)  is  the  preferred  choice 
because  it  can  be  calculated  quickly  using  matrix  algebra.  [22] 
shows  how  a  design  matrix  M  can  be  using  to  quickly  compute 
the  RSS. 

The  algorithm  initializes  by  randomizing  a  solution  which 
satisfies  equations  6  and  7.  The  initial  parameters  ay  and  /3j 
are  calculated  as  well  as  the  objective  score,  RSS0.  In  swap 
k  for  scan  t  choose  i,l  G  {1,...,P}  detections  and  j,m  G 
{1, . . .  ,P}  targets  such  that  y}t]  =  1  and  yftm  =  1.  Switch 
the  detection  association  so  that  t/^+t  =  \  ancj  y}^  =  1. 


Compute  aj,/3j,am,  (3m,  and  RSSk.  If  the  objective  score 
improves,  the  swap  is  kept,  otherwise  it  is  rejected.  The 
algorithm  then  advances  to  the  next  scan  where  the  same 
process  is  repeated,  and  it  terminates  once  it  makes  a  single 
pass  through  all  scans  without  accepting  a  single  switch.  As  we 
will  see  in  the  computational  results  section.  Algorithm  1  runs 
very  efficiently,  providing  high  quality  global  solutions  very 
quickly.  Furthermore,  this  algorithm  can  be  parallelized  by 
running  partitions  of  the  N  starting  points  on  separate  cores, 
leading  to  even  greater  performance  advantages.  A  proposed 
pseudocode  for  the  heuristic  is  provided  below  in  Algorithm  1 . 

Algorithm  1  Randomized  local  search  with  heuristic  swaps 

Input:  X,  P,  T,  M 

Output:  RSS 

Initialization  :  Assign  random  initial  assignments  for  y°tj 
1 :  Calculate  ay ,  /3y  Vj 
2:  Calculate  RSS0 
3:  swapped  true 
4:  k  i —  1 

5:  while  swapped  do 
6:  swapped  -s—  false 

7:  for  t  in  {fi,f2,  do 

8:  Randomly  choose  j,  m  G  {1, . . . ,  P} 

9:  Find  i,l  such  that  ykff}  1  and  y'jff 1  «—  1 

10:  Swap  such  that  ykt:/  1  and  yktm  -f—  1 

li:  Calculate  RSSk ,  ctj,  (3j,  am,  (3m 

12:  if  ( RSSk  >  RSS^1)  then 

13:  yk  g-  yk~1 

14:  else 

15:  swapped  g-  true 

16:  end  if 

17:  end  for 

18:  k  G-  k  +  1 

19:  end  while 

20:  return  RSSk ,  ykt] 

V.  Robust  MIO  Model 

In  this  section  we  treat  the  case  of  detection  ambiguity.  The 
key  difference  is  that  now  the  number  of  targets  is  unknown, 
and  this  becomes  a  third  problem  which  we  wish  to  solve 
in  addition  to  the  data  association  and  trajectory  estimation 
problems  which  remain  once  the  number  of  targets  has  been 
determined.  Since  in  this  case  both  missed  detections  and  false 
alarms  are  present  the  number  of  targets  is  unknown  and  we 
may  no  longer  have  the  same  number  of  detections  at  each 
scan.  Therefore,  we  must  introduce  additional  notation  for 
this  scenario.  We  let  nt  represent  the  number  of  detections 
at  time  t.  We  can  then  identify  the  fewest  and  largest  number 
of  detections  in  a  scan  with  Nn  =  min  m  and  TVi  =  max  n+, 

t  t 

respectively. 

Specifically,  in  this  case  we  replace  Assumption  2  by  the 
following  less  restrictive  assumptions. 

Assumption  3.  (i)  The  sensor  does  not  generate  a  detection 

for  any  target  for  any  time  with  probability  Pd  which  is 
constant  and  independent  between  targets  and  scans. 
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(ii)  At  each  point  in  time  the  sensor  generates  false  alarms 
according  to  a  Poisson  distribution  with  rate  A  fa,  which 
are  located  uniformally  in  the  space. 

(Hi)  The  number  of  true  targets  P  satisfies  N0  <  P  <  N\. 

We  first  show  that  this  problem  can  be  solved  by  dividing 
it  into  a  subset  of  simpler  problems.  We  present  a  MIO 
formulation  that  assumes  a  fixed  number  of  targets.  This 
formulation  allows  us  to  leverage  the  power  of  parallelization 
to  solve  the  problem  by  solving  each  subproblem  separately. 
The  results  can  then  be  gathered  and  compared  to  find  the 
globally  optimal  solution.  For  completeness  we  also  present 
a  formulation  which  solves  the  original  problem  without  the 
need  for  multiple  parallelized  MIOs. 


A.  Fixed  Number  of  Targets  (P) 

If  we  first  assume  that  the  number  of  targets  is  fixed,  we 
can  more  easily  adapt  the  generalized  formulation  presented  in 
Section  III  to  handle  the  addition  of  false  alarms  and  missed 
detections.  This  simple  adaptation  requires  the  introduction  of 
two  additional  variable  types  and  minimal  constraint  changes. 
We  can  then  run  these  formulations  for  each  possible  value 
of  fixed  number  of  targets  P  across  the  range  of  Nq  to 
Ni  and  choose  the  solution  with  the  best  objective  overall. 
Furthermore,  this  is  an  advantageous  strategy  because  each 
independent  experiment  can  be  run  in  parallel. 

1 )  Decision  Variables:  We  first  introduce  new  binary  deci¬ 
sion  variables  Fu  to  indicate  whether  or  not  a  detection  xlt  is 
a  false  alarm. 


Fit 


1,  if  detection  i  at  time  t  is  a  False  Alarm, 
0,  otherwise. 


Similarly,  we  introduce  binary  decision  variables  Mjt  to 
indicate  whether  or  not  an  existing  trajectory  j  has  a  missed 
detection  at  time  t. 


1, 


Mjt  = 


0, 


if  detection  for  trajectory  j 
at  time  t  is  a  Missed  Detection, 
otherwise. 


2)  Constraints:  All  detections  must  either  be  assigned  to  a 
trajectory  j  or  a  false  alarm. 


'y '  Uitj  +  Fa  —  i  vf .  t 
j= i 


(20) 


All  trajectories  j  must  either  be  assigned  a  detection  or  a 
missed  detection. 


E  Vitj  +  Mjt  =  1  Vj,  t 


(21) 


i= 1 


The  sum  of  all  false  alarms  is  TF,  and  similarly  the  sum  of 
all  missed  detections  is  TM. 


nt  T 

i= 1  t= 1 


(22) 


P  T 

EE  Mjt  =  TM  (23) 

o=i  t= 1 

3)  Objective  Function:  We  can  easily  extend  (12)  to  ac¬ 
count  for  false  alarms  and  missed  detections  by  introducing 
penaties  tpo  respectively)  for  each  missed  detection  (false 
alarm,  respectively).  Such  an  objective  form  would  take  the 
form  of: 

p  T 

minimize:  E  E  +^TF  +  fi0TM  (24) 

Zjt,<Xj,Pj,TF,TM  '  j— ' 

0— 1  *=1 

which  can  be  linearized  in  the  same  manner  as  (13). 

4)  Formulation  2: 

P  T 

minimize:  V''  ipjt  +  9aTF  +  fifTM 

*** 

p 

subject  to:  E  Vitj  +  Fu  =  1  Vi,  t 
o= i 
nt 

E  ynj  +  Mp  = 1  vj,  t 

i= 1 
nt  T 

EE^  =  tf 

i=l  t= 1 
P  T 

EEm^  =  ™ 

3=1  t= 1 

■t'it ijitj  F  Mt  ( 1  Vitj )  f.  Zjt  Vi,  t,j 

■t-'it Vitj  Mt(  1  Vitj)  f  Zjt  Vi ,  t.  j 

zjt  -  atj  -  fijt  <  f)jt  Vj,  t 
-  ( zjt  -  otj  -  fijt)  <  ipjt  Vj,  t 
Vitj&{  0,1}  Vi,  t,  j 

OLj  G  K",  P j  G  R"  Vj 

zjteM.n,  Vj,  t 


B.  Number  of  Targets  as  a  Decision  Variable 

In  the  previous  section,  we  assumed  we  knew  the  number 
of  targets.  In  this  section,  the  number  of  targets  is  determined 
via  optimization. 

1 )  Decision  Variables:  Toward  this  goal,  we  introduce  a 
new  binary  decision  variable  Wj  to  indicate  whether  or  not 
trajectory  j  corresponds  to  an  existing  target. 


if  trajectory  j  exists, 
otherwise. 


2)  Constraints:  Most  constraints  remain  similar  to  their 
original  counterparts,  except  now  we  must  account  for  the 
possibility  that  some  trajectories  may  not  exist.  Therefore, 
where  before  we  summed  over  P,  we  will  now  be  summing 
over  Aj .  This  affects  two  constraints. 

All  detections  must  either  be  assigned  to  a  trajectory  j  or 
a  false  alarm.  This  can  be  implemented  similarly  to  (20), 
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except  now  we  sum  over  Ni  because  the  number  of  targets  is 
unknown  but  limited  by  Ni. 

JVi 

e  yuj +fu  =  i  v«,  t  (25) 

7  =  1 


Similarly,  (23)  must  be  adjusted  to  sum  over  the  maximal 
number  of  targets  allowed  Ni  . 

Ni  T 

EE  Mjt  =  TM  (26) 

7  =  1  t=i 


All  existing  trajectories  must  either  be  assigned  a  detection 
or  a  missed  detection. 

nt 

E  W  +  =  W3  v?, 1  (27) 

i=l 


We  restrict  a.?  and  f3j  to  be  zero  if  trajectory  j  does  not 
exist.  This  ensures  only  existing  trajectories  are  penalized  in 
the  objective  function. 

\otj\  +  \Pj\<M0Wj  Vj  (28) 


Since  Nq  <  P  <  Aj ,  we  can  set  Wj  =  1  for  all 
j  =  1, ,  N0,  which  leaves  us  with  only  —  N0  additional 
binary  variables.  We  simply  need  the  additional  constraint 


u>n0+i  >  >  wNl  (29) 


which  guarantees  a  unique  w  solution.  Furthermore,  we  can 
replace  (27)  with  the  following  two  constraints: 

nt 

E  Vit 7  +  M3 1  =  1  Vj  =  1, ...,  No,  t  (30) 

i=  1 
nt 

E  yuj  +  Mjt  =  wj  Vj  =  JV0, iV1;f  (31) 

i=l 


3)  Formulation  3:  Incorporating  these  additional  variables 
and  constraints,  we  arrive  at  the  following  complete  alternative 
formulation. 


Ni  T 

minimize:  |z,-t  —  a7-  —  /3,-fl  +  0OT'-F1  +  4>({FM 

Vitj  ,<Xj  , Fit, Mjt  t=\ 

N\ 

subject  to:  E  Vitj  +  Fu  =  1  Vi,  t 
7  =  1 
nt 

E  Vitj  +  Mjt  =  1  Vj  =  1, ...,  No,  t 

i= 1 
nt 

E  Wi  +  M7*  =  Vj  =  iV0, ...,  Ni,t 

i= 1 
nt  T 

YY.f«  =  tf 

i=  1  £=1 

A/i  T 

EE%=™ 

7  =  1  t=l 

W7V0  +  1  >  •••  >  WATi 

I  ay- 1  +  l^jl  <  Vj 

xuUitj  V  3/|  ( 1  Vitj')  V  7,' 7 1  Vi,  f ,  j 

3'itlhtj  Mi(l  Vitj)  —  Zjt  Vi,  t,j 

yitj&{ 0,1}  Vi,  f,  j 

aj  G  Kn,  , 8j  G  R",  Wj-  G  Vj 

%Gl",  Vj,  i 

C.  Robust  Extension  to  Algorithm  1 

The  heuristic  for  the  scenario  with  ambiguity  follows  simi¬ 
larly  from  the  heuristic  developed  under  the  scenario  without 
ambiguity.  The  main  difference  is  that  now  the  options  for 
making  switches  must  include  false  alarms  and  missed  detec¬ 
tions.  Therefore,  the  framework  of  the  new  algorithm  is  the 
same  as  for  Algorithm  1,  but  the  new  variant  of  the  heuristic 
randomly  chooses  from  the  following  options: 

1)  Switch  detection  assignments  between  two  existing  tar¬ 
gets. 

2)  Switch  the  detection  assignment  of  an  existing  target  with 
a  false  alarm. 

3)  Switch  the  detection  assignment  of  an  existing  target  with 
a  missed  detection  for  a  different  existing  target. 

4)  Move  the  detection  assignment  of  an  existing  target  to  a 
false  alarm  and  replace  it  with  a  missed  detection. 

5)  Move  a  false  alarm  into  the  location  of  a  missed  detection 
for  an  existing  target. 

We  refer  to  this  robust  extension  to  Algorithm  1  as  Algo¬ 
rithm  2.  Similar  to  Algorithm  1,  this  robust  extension  will 
accept  the  switch/move  if  the  objective  score  improves,  and 
reject  the  switch/move  otherwise.  Algorithm  2  terminates  un¬ 
der  the  same  conditions  as  Algorithm  1.  We  expect  Algorithm 
2  to  run  slightly  slower  due  to  the  increase  in  potential 
combinations  of  solutions. 

VI.  Scenario  Complexity  &  Performance  Metrics 

There  does  not  exist  a  unified  approach  for  measuring 
scenario  complexity  as  stated  by  [1]  nor  does  there  exist  clear 
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measures  of  performance  for  each  of  the  trajectory  estimation 
and  data  association  problems.  In  this  paper,  we  argue  that  the 
data  association  problem  has  a  natural  performance  metric  but 
lacks  a  measure  of  complexity,  while  the  trajectory  estimation 
problem  has  a  natural  measure  of  complexity  but  lacks  a  clear 
performance  metric. 

In  the  case  of  the  data  association  problem,  the  preferred 
performance  metric  often  used  in  practice  is  %  accuracy  i.e., 
the  number  of  correct  detection  assignments  out  of  the  number 
of  possible  correct  assignments.  For  the  case  without  sensor 
ambiguity,  the  number  of  possible  assignments  is  simply  the 
total  number  of  detections,  or  equivalently,  the  number  of 
targets  multiplied  by  the  number  of  scans. 


#  correct  assignments  #  correct  assignments 

Accuracy  = -  =  - — - 

Total  #  of  detections  PT 

(32) 

In  the  case  of  sensor  ambiguity,  however,  the  number  of 
possible  correct  assignments  requires  a  deeper  explanation.  To 
develop  a  better  understanding,  we  consider  our  goal,  which  is 
to  correctly  assign  detections  to  targets  and  identify  both  false 
alarms  and  missed  detections.  With  this  in  mind,  we  define 
the  number  of  possible  correct  assignments  as  the  number  of 
targets  multiplied  by  the  number  of  scans  plus  the  number  of 
false  alarms 


#  correct  assignments 

Accuracy  =  — — - . 

PT  +  #  False  Alarms 


(33) 


Whereas  accuracy  serves  as  a  good  measure  of  performance 
for  data  association,  there  does  not  exist  a  corresponding 
measure  of  complexity  which  comparatively  measures  the 
difficulty  of  the  data  association  problem.  We  argue  that 
cr  alone  is  not  the  best  measure  of  difficulty  for  the  data 
association  problem.  For  example,  a  scenario  with  very  close 
target  trajectories  may  not  actually  be  difficult  to  ascertain  data 
associations  even  for  small  a  values,  and  similarly  with  high 
enough  a  values  even  widely  spaced  targets  could  be  difficult 
to  differentiate.  Therefore,  we  introduce  a  metric  p  to  quantify 
this  complexity.  For  ease  of  notation  in  developing  this  metric 
we  first  define  l),jt  as  the  distance  between  one  true  trajectory 
i  and  another  true  trajectory 


p ,  which  is  the  proportion  of  detection  pairs  that  fall  within  a 
closely  defined  proximity  to  each  other. 


X  X)  °ijt 


t=r i<j 


(35) 


This  metric  has  several  desirable  attributes.  First  and  fore¬ 
most,  it  falls  within  the  range  of  [0, 1],  identical  to  the  range 
of  accuracy,  making  it  easily  comparable.  Secondly,  it  is  easy 
to  understand  and  interpret.  Higher  values  of  p  indicate  easier 
scenarios  because  fewer  targets  are  within  close  proximity  for 
a  shorter  amount  of  time,  and  vice  versa.  Finally,  as  we  have 
defined  it,  p  has  an  inverse  relationship  with  cr,  which  means 
that  it  serves  as  a  connection  between  scenario  generation 
and  performance  measuring  processes.  While  a  can  be  used 
more  naturally  for  scenario  generation,  where  it  is  useful  as  a 
parameter  for  signal  noise,  p  can  be  calculated  after  the  fact 
and  used  to  quantify  the  difficulty  of  the  scenario  as  it  pertains 
to  the  data  association  problem. 

In  the  case  of  the  trajectory  estimation  problem,  the  pre¬ 
ferred  complexity  metric  often  used  in  practice  is  cr.  Increasing 
the  signal  noise  may  often  lead  to  stronger  bias  in  the  trajectory 
estimation,  especially  in  scenarios  with  fewer  scans,  and 
results  in  a  deteriorated  quality  of  the  estimation.  Therefore, 
we  believe  that  a  is  the  correct  metric  for  use  in  measuring 
the  difficulty  of  the  trajectory  estimation  problem. 

However,  establishing  a  performance  metric  for  the  trajec¬ 
tory  estimation  problem  is  necessary.  We  choose  to  implement 
a  metric  which  captures  the  core  goal  of  the  trajectory  esti¬ 
mation  problem,  that  is  to  estimate  a  trajectory  as  close  as 
possible  to  the  true  ground  track. 


6  = 


T  P 


X  X  \\Xjt-Xjt\\ 

t=lj=l 


PT 


(36) 


We  match  the  true  trajectories  to  the  estimated  trajectories 
using  a  one-to-one  assignment  problem  which  can  be  formu¬ 
lated  using  linear  optimization.  Lower  values  of  S  correspond 
to  higher  performance  because  the  distance  between  the  esti¬ 
mated  and  true  ground  trajectories  is  smaller. 

In  Section  ,  we  will  see  how  these  measures  of  complexity 
and  performance  are  useful  in  quantifying  the  strengths  and 
weaknesses  of  our  methods. 


Dijt  =  || alrue  +  firuet  -  afue  +  /3jruet\\  (34) 

Additionally,  we  define  a  variable  Cijt.  that  will  take  the 
value  of  1  if  the  distance  between  trajectory  i  and  trajectory 
j  is  greater  than  some  constant.  We  propose  the  use  of  2 a, 
since  it  is  hard  to  distinguish  detections  which  lie  between 
target  trajectories  closer  that. 

1,  if  D^t  >  2cr, 

0,  otherwise. 

Then  the  difficulty  of  a  scenario  in  the  sphere  of  the  data 
association  problem  is  quantified  by  the  complexity  measure 


VII.  Experimental  Simulations 

There  does  not  exist  among  the  literature  a  clearly  defined 
comprehensive  set  of  simulation  scenarios  as  pointed  out 
by  [1],  However,  this  work  also  noted  that  two  types  of 
scenarios  of  particular  importance  include  crossing  trajectories 
and  parallel  trajectories.  In  agreement,  we  choose  to  develop 
scenarios  of  both  types.  We  evaluated  our  methods  on  two 
separate  experiments,  one  with  detection  ambiguity  and  one 
without. 

Both  experiments  were  implemented  in  the  development 
software  julia  0.4.3  [23]  using  the  optimization  package  JuMP 
[24],  The  implemented  MIO  utilized  the  optimization  solver 
Gurobi  6.5.0[25],  Gurobi  was  limited  to  the  use  of  a  single 
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core  for  the  optimization  processes.  Each  simulation  was  run 
on  a  single  node  of  the  unclassified  TX-Green  cluster  located 
at  Lincoln  Laboratories  [26].  The  cluster  utilizes  DL165  G7 
compute  nodes,  consisting  of  2.2  GHz  compute  nodes,  with 

8  GB  of  RAM  each,  for  a  total  peak  performance  of  77.1 
TLLOPS. 

A.  Experiment  1 

In  order  to  evaluate  scalability  we  test  our  methods  across 
a  range  of  scenarios  with  varying  numbers  of  targets  and 
scans.  In  particular  we  consider:  P  £  {4, 6, 8, 10}  and 
T  £  {4,  6,  8, 10}  seconds.  Scans  are  collected  at  a  rate  of 
1  Hz.  The  cartesian  product  of  P  and  T  creates  16  unique 
scenario  sizes.  We  generate  10  unique  crossing  scenarios  and 
10  unique  parallel  scenarios  of  each  size.  Crossing  scenarios 
have  trajectories  that  intersect  through  time,  while  parallel  sce¬ 
narios  have  trajectories  within  close  proximity  to  one  another 
but  do  not  ever  actually  intersect.  Trajectories  are  restricted 
to  exist  within  a  fixed  positional  window  of  [—10,10].  Lor 
each  scenario,  we  randomly  generate  10  realizations  of  data 
by  perturbing  each  true  position  measurement  by  an  error 
e  ~  7V(0,  cr)  with  cr  e  {0.1,  0.5, 1.0,  2.0,  3.5,  5.0},  where 
a  represents  the  noise  parameter.  The  problem  data  is  then 
generated  by  adding  the  detection  error  to  the  true  position. 

_  true  .  ptruej.  ,  _  soh\ 

Xit  =OLi  +  Pi  t  +  e  (37) 

Scans  Xf  are  simulated  by  randomizing  the  order  of  xit  for 
each  t.  Each  unique  X  generated  is  referred  to  as  a  simulation. 
Lor  each  such  simulation,  we  run  the  heuristic  with  a  range  of 
starting  points  N  £  {100  1,000  10,000},  and  use  each  of 
these  solutions  as  a  warmstart  for  the  MIO.  The  optimization 
process  is  set  to  terminate  after  3T  seconds,  with  solutions 
collected  at  intervals  of  {1,  T,  2T,  3T}  seconds. 

B.  Experiment  2 

The  second  experiment  serves  as  an  extension  of  the  first 
in  order  to  test  the  performance  of  our  algorithms  under 
detection  ambiguity.  We  use  the  same  base  data  generated  from 
Experiment  1,  but  now  we  simulate  missed  detections  and  false 
alarms.  A  detection  is  removed  with  probability,  7,  and  we 
consider  7  £  {0.8,  0.85,  0.9, 0.95}.  Lor  each  scan,  we  generate 
false  alarms  by  a  poisson  distribution  with  parameter.  A,  which 
locations  are  then  randomly  selected  uniformly  within  the 
state  space.  The  false  alarms  are  then  added  to  Xt  and  the 
detection  order  of  Xt  is  randomly  shuffled.  Once  the  data 
has  been  generated,  we  use  the  same  approach  as  Experiment 
1,  testing  our  algorithms  with  identical  values  of  N  and 
capturing  solutions  at  the  same  intervals.  Additionally,  we 
test  the  sensitivity  of  our  methods  across  a  range  of  penalties 

9  £  {TBD}  and  0  £  {TBD}. 

VIII.  Computational  Results 

We  begin  by  discussing  the  relationship  between  p  and  a 
and  discuss  the  how  this  relationship  benefits  both  scenario 
generation  and  measuring  complexity.  Then  we  frame  the 
performance  of  the  basic  heuristic  before  discussing  the  perfor¬ 
mance  of  the  basic  MIO  model  in  both  the  data  association  and 


trajectory  estimation  spheres.  We  follow  this  with  a  discussion 
of  the  robust  MIO  model  evaluated  under  both  spheres. 

A.  Scenario  Generation 

Ligure  1  below  shows  the  relationship  between  a  and  p. 
The  plot  is  broken  down  by  scenario  type  between  crossing 
and  parallel  trajectories.  It  can  be  seen  that  according  to  p, 
the  parallel  method  of  scenario  generation  on  the  average 
creates  easier  scenarios  for  the  data  association  problem.  This 
trend  would  be  expected  because  it  reasons  that  crossing 
scenarios  would  be  more  likely  to  exhibit  detections  within 
close  proximity.  Lrom  this  plot  we  also  see  how  a  small  range 
of  six  values  of  a  corresponds  to  the  full  range  of  p  from  0  to 
1,  meaning  that  we  can  quantify  data  association  performance 
across  a  more  continuous  range. 


Sigma  vs  Rho  for  20  Scenarios  by  Scenario  Type 


o.oo- .  . 

0.1  0.5  1  2  3.5  5  0.1  0.5  1  2  3.5  5  0.1  0.5  1  2  3.5  5  0.1  0.5  1  2  3.5  5 

Sigma 


Lig.  1:  Relationship  between  cr  and  p  summarized  by  scenario 
type. 


B.  Basic  Heuristic  Results 

We  begin  our  discussion  of  the  heuristic  by  evaluating  the 
run  times.  Table  I  summarizes  the  average  run  times  of  the 
heuristic  from  Experiment  1  for  each  value  of  N,  the  number 
of  heuristic  starting  points,  arrange  by  the  number  of  targets 
(P)  and  number  of  scans  ( T ). 

We  see  that  run  times  increase  linearly  with  the  number 
of  starting  points  for  a  fixed  number  of  targets  and  scans. 
Because  of  this  relationship.  Table  I  could  be  used  to  easily 
compute  the  run  time  required  for  a  given  application.  Lor  a 
given  scenario  of  targets  and  scans  and  a  desired  number  of 
starting  points,  the  required  heuristic  run  time  can  be  projected 
using  this  table.  While  this  may  seem  to  grow  quickly  for 
some  scenarios,  it  is  important  to  note  that  the  power  of  the 
heuristic  is  that  it  can  be  leveraged  through  parallelization  by 
running  a  subset  of  the  total  desired  number  of  starting  points 
on  several  processors.  aThe  linear  relationship  works  to  our 
advantage  with  parallelization  because  the  total  mn  time  for  all 
starting  points  can  be  divided  by  the  number  of  processors  to 
be  used  in  parallelization  to  find  the  required  run  time  for  each 
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p 

T 

Number  of  starting  points 

100 

1,000 

10,000 

4 

4 

0.01 

0.10 

0.86 

4 

6 

0.02 

0.24 

2.18 

4 

8 

0.05 

0.45 

4.28 

4 

10 

0.08 

0.76 

7.34 

6 

4 

0.01 

0.15 

1.41 

6 

6 

0.04 

0.39 

3.67 

6 

8 

0.09 

0.81 

7.87 

6 

10 

0.17 

1.56 

15.34 

8 

4 

0.02 

0.19 

1.83 

8 

6 

0.07 

0.57 

5.49 

8 

8 

0.14 

1.24 

12.19 

8 

10 

0.27 

2.53 

25.15 

10 

4 

0.03 

0.25 

2.38 

10 

6 

0.09 

0.80 

7.87 

10 

8 

0.20 

1.84 

18.48 

10 

10 

0.39 

3.73 

37.44 

TABLE  I:  Heuristic  run  times  on  a  single  processor  (in 
seconds). 


subset  of  starting  points.  For  example,  consider  a  scenario  of 
six  targets  over  a  period  of  six  scans.  If  we  desire  to  run  the 
heuristic  with  50,000  starting  points,  the  total  run  time  required 
when  run  in  sequence  is  projected  to  be  approximately  40 
seconds  according  to  table  one.  However,  if  we  parallelize 
this  onto  100  processors,  then  the  real  run  time  is  reduced 
to  0.4  seconds.  Thus,  the  run  time  of  the  heuristic  can  be 
reduced  to  meet  efficiency  needs  subject  only  to  the  limitation 
of  available  processors. 

We  continue  our  evaluation  of  the  basic  heuristic  by  ana¬ 
lyzing  the  performance  of  its  solution  on  the  MIO  objective. 
This  gives  insight  into  wether  or  not  the  use  of  RSS  as  a  proxy 
for  the  MIO  objective  is  effective.  To  evaluate  effectiveness 
on  these  terms,  we  compute  the  corresponding  MIO  objective 
value  of  the  heuristic  solution  and  normalize  it  against  the 
MIO  objective  value  of  the  ideal  solution,  which  refers  to 
the  solution  in  which  the  data  association  problem  is  exactly 
correct  (all  detection  assignments  are  exactly  known).  We  refer 
to  the  resulting  normalized  value  as  percentage  of  the  ideal 
solution’s  MIO  objective  value.  Figure  2  plots  this  normalized 
against  cr. 

We  see  that  increasing  the  number  of  starting  points  im¬ 
proves  the  quality  of  the  heuristic  solution  as  compared  to 
the  ideal  solution’s  MIO  objective  value,  especially  when  the 
number  of  targets  is  small.  However,  this  effect  is  diminished 
as  the  number  of  targets  increases.  In  addition,  we  see  that  for 
larger  numbers  of  targets,  even  the  largest  number  of  starting 
points  does  not  achieve  near  ideal  performance,  suggesting  the 
need  for  a  much  larger  number  of  starting  points.  This  is  not 
considered  to  be  a  problem,  however,  due  to  the  advantages  of 
parallelization  discussed  previously  and  also  due  to  the  power 
of  optimization,  which  we  will  see  later. 

We  also  see  that  for  larger  values  of  cr  the  heuristic  actually 
outperforms  the  ideal  solution’s  MIO  objective  value.  Remem¬ 
ber  that  the  ideal  solution  is  simply  ideal  in  the  sphere  of  data 
association,  while  the  MIO  objective  intends  to  score  both 
the  data  association  and  trajectory  estimation  simultaneously. 
Therefore,  we  draw  the  conclusion  that  achieving  perfect  data 
association  for  large  values  of  sigma  does  not  necessarily 
correspond  to  the  best  solution  to  the  trajectory  estimation 


Fig.  2:  Heuristic  performance  as  a  percentage  of  the  ideal 
solution’s  MIO  objective  value. 


problem.  In  other  words,  as  cr  increases  it  may  be  necessary 
to  tradeoff  correct  data  associations  in  order  to  improve  the 
trajectory  estimation.  We  believe  the  results  of  the  heuristic 
could  be  explained  by  this  effect. 

Next,  we  evaluate  the  performance  of  the  basic  heuristic  on 
the  data  association  problem  as  a  function  of  the  number  of 
starting  points.  To  this  end,  we  relationship  between  accuracy 
and  N.  Figure  3  plots  the  mean  accuracy  of  each  of  the  three 
starting  points  from  Experiment  1  against  p. 


Mean  Heuristic  Accuracy  vs  Rho  by  Number  of  Starting  Points  (N) 


T:4  T:  6  T:8  T:  10 
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Fig.  3:  Accuracy  of  basic  heuristic  by  number  of  heuristic 
starting  points. 

First  of  all,  we  see  that  the  heuristic  finds  good  solutions 
to  the  data  association  problem,  especially  for  scenarios  with 
fewer  targets,  but  performance  degrades  as  the  number  of 
targets  increases  which  is  expected.  Again  it  is  seen  that 
increasing  the  number  of  starting  points  results  in  minor 
improvements,  and  this  improvement  is  greatest  for  scenarios 
with  fewer  targets.  We  see  a  similar  effect  as  p  increases. 
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We  conclude  that  even  small  values  of  N  produce  moderately 
good  solutions  as  measured  by  accuracy. 

Overall  we  conclude  that  there  is  not  a  significant  difference 
in  heuristic  performance  for  the  range  of  N  values  that  we 
explored.  Therefore  for  simplification  as  we  move  forward  in 
our  analysis,  we  will  restrict  our  discussions  of  the  heuristic 
to  JV  =  1, 000. 


C.  Basic  MIO  Results 

Next,  we  evaluate  the  accuracy  of  the  MIO.  Figure  4  plots 
the  mean  accuracy  of  the  MIO,  initialized  by  the  N  =  1 , 000 
heuristic  solutions,  after  1,T,  and  2T  seconds  against  p.  We 
have  excluded  the  data  for  the  MIO  after  3T  seconds  for  the 
sake  of  clarity  as  it  showed  little  to  no  improvement  over  the 
MIO  after  2T  solution.  For  comparison,  we  have  included  the 
heuristic  (for  N  =  1,000  only)  in  addition  to  a  randomized 
solution,  one  in  which  we  randomly  assigned  detections  to 
targets.  Note  that  in  this  case,  that  the  ideal  solution,  one 
in  which  the  associations  are  exactly  correct,  achieves  an 
accuracy  of  1.0  in  all  cases  so  it  is  not  shown  explicitly. 


Fig.  4:  Accuracy  of  MIO  compared  against  the  heuristic  and 
a  randomized  solution. 


For  scenarios  with  fewer  numbers  of  targets,  the  MIO 
solutions  were  actually  proven  to  be  the  optimal  solution. 
Therefore,  for  smaller  scenarios  with  few  targets,  we  see  that 
the  heuristic  achieves  optimal  and  near  optimal  solutions.  We 
also  see  that  the  easier  the  scenario,  the  more  improvement  the 
MIO  has  over  the  heuristic,  while  in  more  difficult  scenarios 
the  effect  is  diminished.  Furthermore,  it  can  be  seen  that  in 
nearly  all  scenarios,  the  MIO  achieves  its  best  or  near  best 
solutions  after  T  or  fewer  seconds,  suggesting  the  usefulness 
of  the  MIO  as  an  online  algorithm  with  a  sliding  window. 

Next,  we  evaluate  the  performance  of  the  basic  heuristic  and 
MIO  through  the  lens  of  trajectory  estimation.  As  discussed 
previously,  we  are  interested  in  comparing  <5,  our  proxy  for 
ground  track  error,  against  cr,  our  measure  of  difficulty  for 
trajectory  estimation,  in  order  to  analyze  performance  of  in 
the  sphere  of  estimation.  Figure  5  plots  cr  against  6  for  each 


of  the  solution  types.  In  addition  to  the  random  solution  shown 
on  the  previous  plot,  we  also  add  a  comparison  to  the  ideal 
solution,  as  previously  defined. 


SolutionJType 

-  Random 

-  Heuristic 

-  MIO_1_sec 

-  MIO_T_sec 

-  MIO_2T_sec 
Ideal 


Fig.  5:  Trajectory  estimation  performance 


Remember  that  lower  values  of  delta  correspond  to  trajec¬ 
tory  estimations  that  are  closer  to  that  of  the  true  ground  track. 
We  see  that  the  performance  of  the  heuristic  converges  to  that 
of  the  MIO  for  scenarios  with  few  targets,  as  well  as  for  large 
values  of  a.  Additionally,  we  see  that  as  the  number  of  targets 
increases  we  begin  to  see  stronger  improvements  by  the  MIO 
over  the  heuristic.  Interestingly,  we  see  that  for  the  scenarios 
with  the  largest  number  of  targets  and  scans,  the  MIO  after 
one  second  is  not  much  better  than  the  heuristic.  While  the 
MIO  after  T  seconds  provides  significant  improvement  over 
that  of  the  heuristic  and  MIO  after  1  second,  there  is  little 
further  improvement  in  running  the  MIO  for  2T  seconds. 

Again,  we  see  that  in  scenarios  with  only  for  scans  (T  =  4) 
we  see  that  for  larger  values  of  cr  the  heuristic  and/or  MIO 
sometimes  outperforms  the  ideal.  This  is  likely  a  result  of 
limited  data  and  increases  uncertainty  under  high  noise.  As 
the  number  of  scans  approaches  infinity,  the  ideal  solution,  or 
perfect  data  associations,  leads  to  trajectory  estimates  that  are 
closer  and  closer  to  the  true  ground  track.  Put  differently,  as 
more  and  more  data  is  known,  it  becomes  easier  to  estimate 
the  trajectories  even  in  the  event  of  large  noise,  and  so 
the  trajectory  estimates  that  result  from  the  ideal  solution 
converges  to  the  true  ground  track. 


D.  Robust  Heuristic  &  MIO  Results 
Results  to  be  included  in  next  draft. 

IX.  Conclusion  and  Future  Work 

We  presented  a  multi-target  tracking  approach  which  jointly 
solves  the  problems  of  data  association  and  trajectory  esti¬ 
mation.  We  accomplish  this  without  the  need  of  a  trajectory 
bank  nor  the  a  prior  computation  of  trajectory  hypothesis.  We 
demonstrated  that  the  proposed  method  outperforms  for  linear 
trajectories.... 
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1)  Online  algorithm  with  sliding  window  2)  More  complex 
penalties 
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